Mon. Not. R. Astron. Soc. 000, [TH23] (2010) Printed 15 September 2010 (MN JAT^ style file v2.2) 

The Dance of Heating and Cooling in Galaxy Clusters: 
3D Simulations of Self-Regulated AGN Outflows 
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ABSTRACT 

It is now widely accepted that heating processes play a fundamental role in galaxy 
clusters, struggling in an intricate but fascinating 'dance' with its antagonist, radiative 
cooling. Last generation observations, especially X-ray, are giving us tiny hints about 
the notes of this endless ballet. Cavities, shocks, turbulence and wide absorption-lines 
indicate the central active nucleus is injecting huge amount of energy in the intraclus- 
ter medium. However, which is the real dominant engine of self-regulated heating? 
One of the model we propose are massive subrelativistic outflows, probably generated 
by a wind disc or just the result of the entrainment on kpc scale by the fast radio 
jet. Using a modified version of AMR code FLASH 3.2, we explored several feedback 
mechanisms which self-regulate the mechanical power. Two are the best schemes that 
answer our primary question, id est quenching cooling flow and at the same time pre- 
serving a cool core appearance for a long term evolution (7 Gyr): one more explosive 
(with efficiencies '^ 5 x 10~^ — 10~^), triggered by central cooled gas, and the other 
gentler, ignited by hot gas Bondi accretion (with e = 0.1). These three-dimensional 
simulations show that the total energy injected is not the key aspect, but the results 
strongly depend on how energy is given to the ICM. We follow the dynamics of best 
models (temperature, density, SB maps and profiles) and produce many observable 
predictions: buoyant bubbles, ripples, turbulence, iron abundance maps and hydro- 
static equilibrium deviation. We present a deep discussion of merits and flaws of all 
our models, with a critical eye towards observational concordance. 
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1 INTRODUCTION 

A fundamental gap in our understanding of the forma- 
tion and evolution of galaxies and galaxy clusters con- 
cerns the thermal ev olution of the baryonic co mponent 
of these systems (e.g . iMcNamara fc NulsenI |2007| [MN07]; 
ICattaneo et al.M2009l ). Massive dark matter halos contain 
large amount of hot gas, shining in the X-ray band. The ob- 
served radiative losses, if not compensated for by some kind 
of heating, would imply gas cooling rates ranging from ~ 1 
M0 yr~^, for massiv e elliptical galaxies, to hundre ds M© 
yr"^ for rich clusters (|Fabianlll994l : IPeres et al.ll iggsl ) . How- 
ever, since the first XMM-RGS observations it has been clear 
that the (radiative) cooling rate in clusters and galaxies is 
reduced by at least one order of magnitude with respect to 
the simple expectation (jPeterson et al.ll200ll . 120031 : IXu et al.l 
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|2002| : IPeterson fc Fabianll2006l . and references thereii 
is the so-called 'cooling flow problem'. 

Active galactic nuclei (AGN) can easily provide enough 
energy to the gas to offset the energy lost by radi- 
ation and high resolution X-ray images show indeed 
clear evidenc e of AGN-gas interaction in many clusters 
and galaxies (jBoehringer et al. 19931: iBlariton et al.l I2OOII : 



iFinoguenov fcJonej|200ll : I Jones et "811120021 : MN07 and ref- 
erences therein) . The fairly common presence of X-ray cavi- 
ties, often coincident with lobes of radio emission connected 
to the core of the central galaxy by a radio jet, indicates 
that AGN inject energy in the intracluster medium (ICM) 
in kinetic form (outflows) and as relativistic particles, al- 
though the quantitative significance of the latte r is difficult 
to estimate (|Dunn et al.ll2005l : lGitti et al.ll2009l '). 

iRaffertv et all (|2006l ) and lRaffertv et al.l (|2008l ) showed 
that the AGN power associated to the cavity formation is 
of the same order as the core X-ray luminosity, for a sample 
of 33 clusters and groups. Although this energetic balance is 
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only a necessary but non sufficient requirement for a heating 
scenario to be successful, it strongly suggests that the heat- 
ing process manifests itself generating bubbles in the ICM. 
Cavities can be easily created by 'directional' input of en- 
ergy, such as jets or coUimated outflows, making spherically 
symmetric form of heating less appealing as major players 
in solving the cooling flow problem. 

In order to prevent significant gas cooling, the feed- 
back process must be activated with a frequency not greatly 
different from 1/fcooi, where tcooi is the central cooling 
time, often of the order of few x 10^ y r for clusters (e.g., 
ISanderson et al.ll2006l : lMittal et al.ll2009l ') and even lower for 
elliptical galaxies or groups (e.g. tcooi ~ 1-5 x 10^ yr for NGC 
4636, see lBaldiet"alIl2009l ). 

Moreover, the feedback heating must prese rve the cool 
core appearance of the majority of the clusters (| Peres et al.l 
ll998l : lMittal et al.ll2009l ). In fact, it has been shown that con- 
centrated heating while very efficient in stopping the cool- 
ing process, often generates ne gative temperature gradients , 
contrary to the observatio ns iJBrighenti fc Mathewa l2002l . 
I2OO3I : iMathews et al.ll2006l '). Another indication that AGN 
cannot deposit most of its energy in the very central region 
is the com mon survival of galacti c scale cool cores in cluster 
ellipticals i|Sun et al.ll2005l . 120071 ): a spatially concentrated 
heatin g would easily erase these fragile l ow temperature re- 
gions JBrighenti fc Mathewd[200i . l2003l ). 

Motivated by these considerations, we investigate here 
the long t erm effect of kinetic feedback on the ICM (see 
MN07 and lFabianll2009l '). We assume that AGN outbursts 
generate collimated, subrelativistic outflows on kpc scale. 
There is widespread observational evidence for winds orig- 
inating in galactic nuclei. High redshift radio galaxies host 
galactic scale, bipolar outflows of ionized gas with veloc- 
ity ~ 1000 km s~^, likely trig gered by the i nteraction of 
the radio jet with the ISM (Nesvadba et al.1 ^2008). Mor- 
ganti et al. (2005, 2007) report the detection of fast mas- 
sive neutral outflows, using 21-cm blueshifted absorption 
lines against strong radio continuum. They occur at kpc 
distance from the nucleus with rates of tens M© yr~^. Op- 
tical, UV and X-ray observations of highly ionized gas also 
point toward fast outflows (thousands km s~^) driven by 
entr ainment, even if t h e mass rates are relatively mod- 
est dGeorge et al.1 1 19981: ICrenshaw et al.lll999l: iKrisd l2003l: 



Risaliti e t al. 2005:'McKernan e t al.ll2007J : [Pounds fcReevesI 
2OO9I : see [Crenshaw et aL 2003 for a review). It is reasonable 
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to conclude that most of the AGN exhibit outflows, albeit 
the physical parameters of the winds are still uncertain. The 
geometry of the flow is also unclear, both polar winds, per- 
haps caused by entrainment of the IC M in the re lativistic 
radio jet, or equatorial disc winds f see IProgal 120071 and ref- 
erences therein) being possible. 

The key question we want to address in this work is the 
following: are outflows from the central AGN able to prevent 
the ICM from cooling and at the same time preserve the cool 
core appearance? 

In recent years a considerable amount of research 
has been devoted to the understanding of the effect of 
jets on the ICM. Most works inves tigated the transient 
flow r esulting from a bipolar outflow ( Revnolds et al.l 2001 , 
2002; Basson & Alexander '2003; Ruszkowski et al. 2004; 
Omnia et al. 2004 ; Zanni et al. 2005 ; Brueggen et al. 2007, ; 



[Sternberg erai]|2007l : [Sternberg fc Sokeill2009l ). but did not 
fully confronted the question above. 

The long term influence of AGN outflows on cluster 
cooling flows was studied by iBrighenti fc Mathewa (2006, 
hereafter BM06). They used 2D simulations and a mechan- 
ical feedback scheme, self-regulating injection time (but not 
velocity and power). They calculated the gas cooling rate 
and the azimuthally averaged density and temperature pro- 
flies. It was found that some intermittent bipolar outflow 
models, with velocity in the range 5 x 10'^ — 10* km s~^, 
could shut down gas cooling for many Gyr, while preserving 
the cool core appearance of the cluster. 



ICattaneo fc Tevssien (|2007|) performed 3D calculations 
of AGN feedback in a poor cluster, about ^ 10 time less mas- 
sive than the object studied here. They employed a hybrid 
kinetic-thermal feedback which injects energy at a rate pro- 
portional to the Bondi accretion rate to the central black 
hole. With the adopted ICM initial conditions the central 
cooling time is quite long (4 Gyr), and the average cooling 
rate in absence of feedback is ~ 30 M© yr~^. During the 
most powerful AGN outbursts only a small fraction of en- 
ergy is ejected in kinetic form, the velocity of the jets being 
only ~ 1.4 X 10^ km s~^. Their models are successful in 
stop ping gas cooling, but the cool core disappear after a few 
Gyr. [Dubois et al.l (120_1Q!) implemented this Bondi feedback 
in a cosmological context, with similar results. 

The present paper builds on the preliminary results by 
BM06 to explore a larger set of feedback schemes and bound 
the parameter space for successful feedback models. In this 
work we use multigrid, 3D hydrodynamical simulations for 
two reasons. First, the nature of cooling flows is intrinsically 
chaotic and turbulent, therefore it is essential to allow a re- 
alistic description of all instabilities, which in turn influence 
the outflow evolution. Second, we want also to eliminate the 
spurious cooling of gas along the z symmetry axis present 
in 2D axisymmetric calculations (like in BM06). 

As stated above, our primary objective here is to in- 
vestigate if massive, collimated outflows can provide a suit- 
able form of feedback to solve the cooling flow problem. The 
adopted scheme for the generation of jets and their connec- 
tion with the accretion on the black hole is simplified and 
parametrized (as customary in similar studies), and attempt 
to investigate subjects like growth of the black hole or jet 
physics is certainly superficial. We will critically test our 
models, trying to recognize all the positive and problematic 
aspects of this form of feedback, comparing our results with 
a great variety of observational constraints. 

Finally, we mention that other authors simulated the 
feedback from buoyant cavities and the associated shock 
heating, without explicitly including the jets that likely 
gener at e these features in real cl usters (IBrueggen fc Kaisei 



200^; 'Brighenti fc M athewsl l2003l: 

2004; Dalla Vecchia et all l2004l: IBrueggen et al 



Ruszkowski et al 

_^ .__^ 200i 

Sii acki et al.l 120071: IBrueggen fc Scaiinapiecd I2OO9I : 
IMathews fc Brighentil l2008al lbl: iMathewd l2009l : the lat- 
ter three papers considered cavities generated by cosmic 
rays). This and other different kinds of feedback will be 
also explored by us in forthcoming papers. 



© 2010 RAS, MNRAS 000,[TH23] 



AGN Outflows: Regulating the Dance of Heating and Cooling 3 



2 THE COMPUTATIONAL PROCEDURE 

The simulations presented liere were calculated with a n 
highly modified version of FLASH 3.2 (JFrvxell at alJl200d l. 
a 3D adaptive mesh refinement (AMR) public code, which 
solves the hydrodynamic euler equations through a split 
Piecewise-Parabolic Method (PPM) solver, particularly ap- 
propriate to describe shock fronts. It uses the Message- 
Passing Interface (MPI) library to achieve portability and 
efficient scalability on a variety of different parallel HPC sys- 
tems. The simulations were run on 128 processors of IBM 
P575 Power 6 (SP6) at CINECA supercomputing centre. 

We added several source and sink terms to the usual 
hydro-equations solved by FLASH, in conservative form: 

-^ -H V • {pv) = ap, - q-^ + Si.jct , (1) 

CC (•cool 

^+^-{pv®v) + VP = pgoM + S2.jct , (2) 



-J^ + V • [{pe + P)v]^ pv- gnj^,i + ap. 



V 

eo + - 



P=(7-l)p(e-y 



neniA(T, Z) + Sajct , (3) 

(4) 



where p is the gas density, v the velocity, e the specific total 
energy (internal and kinetic) , P the pressure, g^M the grav- 
ity of dark matter, and 7 = 5/3 the adiabatic index. The 
temperature is computed from P and p using state Eq. Q, 
with an atomic weight fi ~ 0.62, appropriate for a totally 
ionized plasma with 25% He in mass. 

Outflow source terms S'i,2,3,jGt (dependent on injected 
density, momentum and mechanical energy, respectively), 
along with the spatial distribution of the source region, will 
be explained for every type of feedback in Section 2.2. Note 
that injection will be done directly into the domain (without 
a mass inflow. Si — 0) or through boundary condition at 
z = {Si> 0). 

Radiative cooling is treated in the code like a source 
term in Eq. ([3}: the ICM loses energy at a volume rate 
ncniA{T, Z), where rio and ni are the number density of 
electrons and ions, and A is th e cooUng function for a tem - 
perature T and metallicity Z ({Sutherland fc Dopitalll993h . 
In the following we a ssume Z = 0.3 Zp, a typical metallic- 
ity for these systems (jTamura et al.ll2004l '). ignoring central 
negative gradients. The minimum temperature allowed is 
10* K. 

We also consider SNIa and stellar winds heating in the 
central elliptical galaxy (p* is the de Vaucouleurs stellar den- 
sity profllc), although their effect is minor in massive clus- 
te rs. The y arc implemented following the parametrization 
of iBrighenti fc Mathewd (|2002h . with a rate (dominated by 
stellar mass loss) a = a, -I- qsn ~ 4.7 x 10~^°{t/tn)~^'^ 
s~^, where tn = 13.7 Gyr is the present time, and with a 
speciflc injection energy £0 dependent on SNIa and stellar 
winds temperature. 

We also add a mass dropout term (see 
IBrighenti fc Mathewsll2002l ) to avoid a clutter of zones flUed 
with cold gas, which jet events could spread at larger radii 
or, without feedback, accumulate in the nucleus. In fact, 
without this term the simulation will generate clouds of very 



cold gas, whose physics (collapse, star formation) cannot be 
described by the hydrocode. Thus, in order to prevent this 
feature, we include a mass sink term —q{T)p/tcooi in Eq. lU 
and drop out the cold gas (at constant pressure). Its effect 
is simply to remove the cold gas from the grid, without 
affecting the hotter flow o r the calculated cooling rate 
IBrighenti fc Mathewsll2002l ). The dimensionless coefficient 
q{T) is defined as 2exp( — (T/Tq)^) and dropout becomes 
significant when T < Tq = 5 x 10^ K. The cooling time is 
assumed fcooi = 5P/2neniA. Note that the total mass of 
cooled gas does not depend on the presence of the dro pout 
term or its functional form IjBrighenti fc Mathewsll2000l ). 

We calculate the ffow evolution for 7 Gyr (large cluster 
formed relatively recently), as opposed to most of the works 
cited in Section 1. As shown in Section 3.1 several feedback 
heating schemes are only able to delay excessive gas cooling 
for 1-2 Gyr before failing, and therefore we stress the need to 
investigate the long term (several Gyr) behaviour of heated 
ffows. 



2.1 The cluster model and initial conditions 

As in BM06, we adopt the well observed cluster Abell 
1795 (JTamura et al.ll200ll : JEttori et al.ll2002l ) as a template 
for our models. Be ing A 1795 a rather typical, relaxed 
(|Buote fc Tsailll996l ). cool core massive (Mvir ~ lO^'^ M©) 
cluster, all the results we present should be relevant for 
any object in thi s category. The sh ort central cooling time, 
^cooi ~ 4 X 10* yr IjEttori et al.ll2002l ) assures that this cluster 
should host a strong cooling flow in absence of appropriate 
feedback. However, recent X-ray observations place an up- 
per limit to the radiative cooling rate Mcoo i ^ 30 M0 yr~^ 
(jPeterson et al.l 120031 : iBregman et al.ll2006l ). which is a key 
constraint for the success of a model. 

We start our calculations with the hot gas in spheri- 
cal hydrostatic equilibrium in the potential well of the dark 
matter halo: 

dP , . , .d(l>-DM . . , , 

_(,).-p(,)__(,), (5) 

where r is the spheric al radius. The dark m atter halo follows 
a NFW distribution IjNavarro et al.lll996r ). with virial mass 
10^^ Mq, and thus a potential given by: 



4>DM{r) = 



GIVU 



ln(l -I- r/ra 



'■s/(cio2) r/rs 



(6) 



where the concentration at overdensity 102 is C102 ~ 
rvir/j's — 6.6 with virial radius ~ 2.6 Mpc, and /(C102) = 
ln(l-|-cio2) — cio2/(l-l-cio2). We have adopted a ACDM cos- 
mological universe with Om ~ 0.27, flA ~ 0.73 and Ho = 71 
km s"^ Mpc"\ 

Combining Eq. ((Sjl and an observed T{r) fit (see third 
panel of Fig. 1, dotted line) with Eq. ([5]), we recover the den- 
sity radial profile p{r), assuming an ideal gas P = k^^pT/nnip 
and a gas fraction of 0.15 at virial radius.. 

In most models we ignore the contribution to the grav- 
itational potential due to the central galaxy. This is im- 
portant only in the inner few tens kpc and has no signifi- 
cant effect on the cooling rate or the hydrodynamical vari- 
able profiles. In one model we have also included the grav- 
ity fro m the centr al galaxy, modelled as a de Vaucouleurs 
profile (jMellier fc Mathezi:1987 ') with total steUar mass M* 
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Table 1. Parameters and properties of all simulated models. 







e 


W'jct 


^jct 


''jot 


Tjct 


Model 


Feedback 


efficiency 


jet width (kpc) 


jet height (kpc) 


jet velocity (km s~^) 


duration/cycle (Myr) 


CP 


no AGN heating 


- 


- 


- 


- 


- 


Al 


AM,„„i 


5 X 10-* 


2.7 


6.8 


(2eAM,„„ic2/Mact)i/" 


self-regulated 


A2 


AMcooi 


10-3 


2.7 


6.8 


(26AM,ooic2/Mact)l/2 


self-regulated 


A2L 


AMcooi 


10-3 


2.7 


17 


(2eAMcooic2/Mact)i/2 


self-regulated 


A3 


AMcooi 


5 X 10-3 


2.7 


6.8 


(2eAM,ooic2/Mact)l/2 


self-regulated 


A3L 


AMcooi 


5 X 10-3 


2.7 


17 


(2eAM,oolc2/Mact)l''2 


self-regulated 


A3S 


Afcool 


5 X 10-3 


2.7 





(2eAMeoolc2/pjetAAt)l/3 


self-regulated 


A4 


AMcooi 


10-2 


2.7 


6.8 


(2eAM,ooic2/M,,t)i/2 


self-regulated 


A4L 


AM,„„i 


10-2 


2.7 


17 


(26AM,„„ic2/Mact)l/2 


self-regulated 


A5 


AM^ooi 


5 X 10-2 


2.7 


6.8 


(2eAM,ooic2/Mact)i''2 


self-regulated 


Bl 


intermittent 


- 


2.7 


6 


8 


10'* 


20/200 


B2 


intermittent 


- 


2.7 


6 


8 


10* 


10/100 


B3 


intermittent 


- 


2.7 


6 


8 


10" 


1/10 


CI 


continuous 


- 


2.7 


6 


8 


2 X 103 


continuous 


C2 


continuous 


- 


2.7 


6 


8 


6 X 103 


continuous 


C3 


continuous 


- 


2.7 


6 


8 


10* 


continuous 


BONDI 


AfB.lOkpc 


10-1 


2.7 


6 


8 


(2eMBAte2/Afact)l''2 


self-regulated 


BONDI2 


A^B,5kpc 


10-1 


2.7 





(2eMBc2/pj,,A)i/3 


self-regulated 



~ 6 X 10^^ M0 and effective radius Vc ~ 8.5 kpc, and veri- 
fied that its influence on the results is inconsequential (see 
Section 3.1.3). 

The computational rectangular box in all of our models 
extends slightly beyond the cluster virial radius r^ii- We sim- 
ulate just half cluster with symmetric boundary condition 
at z = 0, while elsewhere we set prolonged initial conditions 
with only outflow permitted. Despite the AMR capability 
of FLASH, we decided to use a number of concentric fixed 
grids in cartesian coordinates. This ensures a proper resolu- 
tion of the waves and cavities generated in the cluster core 
by the AGN outflows. We use a set of 9 grid levels (with 
basic blocks of 8 x 8 x 4 points), with the zone linear size 
doubling among adjacent levels. The finest, inner grid has a 
resolution of ~ 2.7 kpc and covers a spherical region of 100 
kpc in radius. In general, grids of every level extend radially 
for about 40 cells. The relatively low resolution is due to the 
need of covering large spatial scales (kpc up to Mpc) and at 
the same time integrating the system for several Gyr, using 
moderate computational resources (50,000 CPU hours). 



2.2 Outflow generation 

We adopt a purely mechanical AGN feedback in form of 
nonrelativistic, coUimated outflows (similar to BM06). In 
this paper we show results only for models with cylindrical 
jets, with velocity parallel to the z— axis. We have calculated 
few simulations with conical o utflows (with half-ope ning an- 
gle up to 70 degrees - see also [Sternberg et al.ll2007l ) and we 
have verified that they have a similar impact on the global 
properties of the flow. In fact, the pressure of the ICM col- 
limates the outflows within few tens kpc (see BM06). 

We consider several types of feedback. In feedback 
scheme A an outflow is activated only when gas cools to 
very low temperature within a spherical region r < 10 kpc, 
and drops out from the flow (conceptuall y, this is similar to 
the " cold feedback model" described bv lPizzolato &: Sokeij 
I2OO5I ). Usually most of the gas cools at the very centre of the 



cluster. We assume that at any timestep a fraction e of the 
rest mass energy of the cooled gas, AMcooic^, is injected as 
kinetic energy. Here AMcooi is the gas mass cooled in a given 
timestep of the flnest grid. This energy is given to the hot gas 
located in a small region at the centre of the grid (the 'active 
jet region'), whose size is indicated in Table 1 and contain- 
ing a gas mass Mact (there is no new injected mass in Eq. 
IT])). At every timestep we set the z component of the veloc- 
ity within the active region to «j = (2eAMcooic^/-M"act)"^ , 
since -Bkjet = 0.5Mactiifet = eAMcooiC^- We will see (Section 
3) that the frequency and strength of the feedback events 
strongly depend on the mechanical efficiency e, which has 
typical values lO'* - IO-2. 

We remark again that this scheme to link gas cooling, 
black hole accretion and outfiow generation does not have a 
strong physical basis: it must be taken as a simplified way to 
implement a self-regulated feedback, which triggers heating 
only when it is needed to halt ICM cooling. In considering 
massive slow outflows, we are implicitly assuming that the 
relativistic radio jet entr ains some ICM m ass (Mact)- More- 
over, other authors (e.g. lGiovanninill2004h found that radio 
jets in cluster central galaxies are highly relativistic on pc 
scale, but rapidly decrease to subrelativistic velocities within 
few kpc from the black hole (especially in Fanaroff-Riley I 
sources), because of the interaction with the dense ISM in 
the inner region. 

In feedback method B the outflows are triggered inter- 
mittently, at fixed times (Table 1) and a fixed velocity Ujct 
is given to the gas located in the active region (again 5*1 = 
in Eq. HI). Feedback scheme B is not self-regulated, but the 
AGN outbursts are forced to occur with a frequency (typi- 
cally of the order of 10^ — 10^ yr~'^) which agrees wit h the 
observational estimates (|Sanders fc Fabiar]l2007l . l2008l ). The 
power of any jet event depends only on the mass present in 
the source region and is independent of the accretion rate. 

Next we run a few models where the outflows are con- 
tinuously generated, that is at every time step we keep the 
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2— component of the velocity within the jet active region at 
a given value, listed in Table 1 (scheme C). 

We also adopt a scheme, named BONDI, in which 
the accretion rate is calculated with the Bondi pre- 
scription (see also ICattaneo fc Tevssien 120071 ') : Msondi = 
47r(GMBH)^Po/Cs05 where po is the volume-weighted hot gas 
density calculated within r <, 5 (or 10) kpc, while Cso is 
the mass- weighted sound speed in the same region. Needless 
to say, the Bondi radius, tb ~ 50 pc, is far smaller than 
our resolution limit, so we refrain to attach a strict physical 
meaning to Mb. In this sense the high mechanical efficien- 
cies (0.1) used for Bondi models are due to the fact that the 
accretion should be a factor 10-100 larger (because of higher 
inner po and lower Cso)- This feedback is fundamentally dif- 
ferent from feedback A because it produces a self-regulated, 
quasi-continuous, low power AGN activity, in contrast to the 
few violent AGN outbursts characteristic of the best models 
adopting scheme A. 

Finally, in BONDI2 and ASS the outflow is not injected 
as usual in the active region, but with a mass, momentum 
and energy flux through the boundary at z = (thus Si > 
0), with a square area —1.35 ^ x <y 1.35 kpc, —1.35 ^ 
y ^ 1.35 kpc. The velocity of the jet is then calculated as 
Hjet ~ (2eAffC'^/pjetA)^ , where M{ is the accretion rate of 
the feedback scheme, A ~ 7.3 kpc'^ is the area through which 
the jet is injected, and pjot ~ 2 x 10~^® g cm~"^ (~ 10~^ the 
initial central gas density). We fix the temperature of the 
jet to low values (10^ K) in order to keep injected thermal 
energy on negligible levels compared to the kinetic flux. 



3 RESULTS 

We describe in this Section the results for various flows, 
exploring a large set of feedback parameters. In order to 
understand the long term behaviour of the models we have 
evolved them for 7 Gyr. The numerical resolution adopted 
does not allow a deep study of the disturbances generated 
by the outflows (such as cavities or shocks). On the contrary, 
we believe these models are appropriate to investigate the 
global properties of the flows, such as the cooling rate and 
the azimuthally averaged density and temperature profiles. 
The suite of simulations described here is used to thoroughly 
explore the outflows parameter space in order to bound the 
region of successful models. 



3. 0. 1 Pure cooling flow 

As a reference flow we first calculated a pure cooling flow 
(CF), where no AGN feedback was used, shown in Fig. 1. 

As expected, both density and (mass weighted) temper- 
ature profiles steepen in the central region, an effect caused 
by radiative losses and the consequent subsonic gas infiow. 
In ~ 1 Gyr the calculated profiles disagree with the obser- 
vations of Abell 1795, although not by a great extent. It is 
interesting, however, that the logarithmic slope in the core 
region of the model temperature profiles at late times is very 
close to 0.4, which Sanderson et al. (2006) found to be typical 
in the central region of cool core clusters. The similarity be- 
tween the temperature distribution in our pure cooling flow 
run and real clusters, where heating is currently preventing 
gas cooling, put severe constraints on the feedback process: 




Figure 1. Evolution of model CF (no AGN feedback). In the top 
panel is shown the gas cooling rate versus time. The middle and 
bottom panels show the temporal evolution of the gas (electron) 
number density and mass weighted temperature profiles, respec- 
tively. The profiles are displayed at 15 different times, as indicated 
in the lowermost panel. Observational data of A 1795 are shown 
with fille d circles llTamura et al.ll200ll . XMM-Newton) and filled 
triangles llEttori et al.ll2003 . Chandra). 



it must not greatly perturb the temperature profile shaped 
by r adiative cooling. This is a deman ding requirement (see 
also iBrighenti fc MathewdlJOO^ . l2003l l . 

After few Gyr the fio w reaches an approximate steady 
state (see also lEttori fc Brighenti 2008 for a quantitative de- 
scription of the temporal change in the observable profiles). 
The bolometric X-ray luminosity slowly increases with time, 
from L^ - 1.5 x 10""^ erg s"^ at t = 0, up to ~ 2.3 x lO'*'^ 
erg s~^ at t = 7 Gyr. The growth of the gas density in the 
cluster core is responsible for the increase of Lx ■ 

It is interesting to investigate the global energetic bud- 
get. The internal energy within r^ir decreases by ~ 5 x 10^^ 
erg (that is, of about 1 %), while the potential energy drops 
by ~ 4 X 10®^ erg, considering both the hot gas remaining in 
the grid and the cooled gas at the centre of the cluster. The 
kinetic energy is ~ 2 x 10^* erg and is therefore negligible. 
Thus, energy is radiated away (i?rad ~ 4.5 x 10^^ erg in 7 
Gyr) mainly at the expense of the potential energy of the 
ICM. 

The gas cools in the very center and is removed from the 
computational grid by the dropout term described in Section 
2. The cooling rate increases with time, reaching ~ 250 M© 
yr~^ at the end of the calculation, a blatant discrepancy with 
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the o bservational results (|Tamura et al.1 I2OOII : lEttori et alj 
12002 1). This is of course the so-called cooling flow problem, 
the focus of this work. All the gas in this simulation cools 
in the very centre of the cluster. This findings are in excel- 
lent agreement with BM06, showing that these results do 
not depend on the hydro-code, the coordinate system, the 
symmetry or the numerical resolution adopted. 

In the following we describe the models with jet feed- 
back. 



3.1 Feedback A (AMcooi regulated) 

3.1.1 Model Al, e = 5 X 10"* 

We start illustrating results for feedback A models, where 
the mechanical energy of the jet is linked to the gas mass 
cooled (within r = 10 kpc). In Fig. 2 we show the relevant 
properties for model Al, with efficiency e = 5 x 10~*. The 
width and length of the active region are 2.7 and 6.8 kpc (1 
and 2.5 grid points), respectively. 

The cooling rate, although reduced with respect to the 
CF model, is still too high (Mcooi ~ 75 Mq yr~^ at the end 
of the simulation). The azimuthally averaged, mass weighted 
temperature and density profiles are similar to those of 
model CF, with only a little temporal variation: the weak 
outflows perturb only slightly the temperature profiles. 

In the right column of Fig. 2 are shown the physical 
characteristics of all jet events. In these plots the quantities 
refer to the half-space 2; J5 considered in our simulations. 
In this model the outfiows are activated frequently, because 
their relatively low mechanical power cannot prevent the 
cooling for a long time. 

In the upper panel of the right column is plotted the 
cumulative (mechanical) energy injected by the jets. At f = 
7 Gyr, E^ct ~ 1.5 x 10*"^ erg has been injected in the ICM (in 
the z ^ Q space). Jets become more frequent at late times, 
because of the slow secular decrease of the central cooling 
time, a result of the slight predominance of radiative losses 
over heating. At the end of the simulation ~ 1.6 x 10^^ M© 
have cooled and dropped out of the hot phase. Clearly, if all 
the cooled gas were accreted on the central black hole, as we 
have assumed in our simple feedback scheme, the final black 
hole mass would result far in excess to that of real black 
holes. Of course, we could formally avoid the problem of the 
excessive black hole mass by assuming that only a fraction of 
the cooled gas actually accretes on it with a higher heating 
efficiency. For instance, this model would be identical if we 
assume that only 1% of the cooled gas is accreted by the 
central black hole and the efficiency is increased to 5 x 10~^. 
This degeneracy means that our model is too simple to allow 
a proper investigation of the black hole growth. 

In the remaining panels are displayed the power, veloc- 
ity and mass cooled (within 10 kpc) during every jet event. 
Typically, the gas in an outfiow is ejected with Hjct ~ 10* 
km s-i and power Pjct = 0.5Mactuiit/At ~ 10*'' - lO*'^ erg 
s~^. AMcoid, shown in the bottom right panel of Fig. 2, is 
the mass cooled within 10 kpc during a given timestep, has 
typical values 10* - few 10* Mq. We note that the Edding- 
ton luminosity, Lfidd ~ 1-5 x 1O^*(Mbh/M0) ~ 1.5 x lO*'^ 
erg s~^ for a 10^ M0 black hole, is close to the mechanical 
power of a typical outfiow. 

In summary, model Al (e = 5 x 10~*) resembles a pure 



cooling fiow model, with a cooling rate still too large and 
must therefore be rejected. The next logical step is to in- 
crease the efficiency e in order to reduce the mass of the 
cooled gas and check if more powerful outfiows perturb the 
variable profiles in an acceptable way. 



3.1.2 Model A2, e = 10"^ 

Model A2, with e = 10""^, is not shown here and we limit to 
a brief description of the results. It has good temperature 
profiles, peaked density profile and a cooling rate of ~ 30 
M0 yr~^ at t — 7 Gyr, which is only marginally accept- 
able. The cooling rate for t <, 2 Gyr, however, is < 10 Mq 
yr~^. Evidently, a low efficiency feedback is able to suffocate 
the cooling fiow for several Gyr. Only at late times (f <; 6 
Gyr) the cooling rate becomes too high. This result empha- 
sizes the importance of calculating a model for many Gyr 
to check the long term thermal evolution. With respect to 
model Al, the jet events are more separated in time, espe- 
cially at early times, consistently with the low cooling rate 
at that epoch. The total amount of the energy transferred 
to the ICM is similar to that for model Al, an indication of 
the self-regulation of the feedback process. 



3.1.3 Model A3, 



5 X 10" 



The increase of the efficiency to e = 5 x 10"'^ generates a 
quite successful model (A3) , which we discuss in more length 
(see also Section 4 for an analysis of the fiow dynamics). The 
cooling rate (Fig. [3} is very low at any time (Mcooi ~ 5 M© 
yr~^), a value fully co mpatible with the current observations 
(A/cooi < 30 Mq vr~*: iPeterson et al.|[200^ : iBregman et al.1 
2006). 

The azimuthally averaged density and temperature pro- 
files for t ^ 1 — 2 Gyr are always in good agreement with 
those observed in A 1795. Both the ICM density and tem- 
perature vary somewhat with time, following the outfiow 
cycles as expected, but the cluster always keeps the status 
of 'cool core cluster', even though gas is essentially not cool- 
ing. The shock waves generated by the outfiows (see Section 
4) are not strong enough to significantly heat the gas and 
perturb the global, azimuthally averaged temperature pos- 
itive gradient, although small amplitude ripples are seen in 
the proflle, corresponding to weak shocks associated to the 
jet propagation. These waves are visible up to a distance of 
~ 400 kpc. 

In this model the AGN feedback is triggered less fre- 
quently, and only about 50 jet events occur. The duty cycle is 
~ 6% (of total time). The total energy injected, ~ 1.25x10'^^ 
erg is again of the same order as in the previous simulations, 
and the average power and velocity of the single outfiow is 
therefore larger. This energy must be compared with the to- 
tal energy radiated away, ~ 1.1 x lO®'^ erg (again calculated 
in the half-space z ^ 0). Not surprisingly the two energies 
are of the same order. In fact, if the gas density distribution 
is similar to that of the standard cooling flow model, then 
the energy radiated away is also similar. So the energetic 
balance to stop the cooling requires that AGN provides an 
energy r^ l(f^ erg. The power of the outflows often exceeds 
the Eddington luminosity. While the latter is not strictly 
relevant in this context, the unpalatable large power might 
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indicate that this feedback fails to simulate the real accre- 
tion process at work. 

It is instructive to examine the evolution of the ener- 
getics in the cluster core region (r ^ 50 kpc), the most 
perturbed by the AGN feedback. After every AGN outburst 
the core kinetic energy increases on average by ~ 8 x lO'''^ 
erg (with the most powerful outflows depositing ^ 2 x 10^^ 
erg), which is comparable to the thermal energy content of 
the core region (-Eth,corG ~ 1-5 x 10*'^ erg). However, the ki- 
netic energy generated by an outflow is dissipated in only 
~ few x 10^ yr. Thus, for most of the time the thermal 
energy dominates over the kinetic energy. Following the dis- 
sipation of the kinetic energy, the thermal energy also rises, 
because of the shock heating. Finally, a large fraction of the 
thermal energy gained by the AGN outburst is transformed 
in potential energy, as a consequence of the quasi-adiabatic 
cooling due to the expansion of the ICM. 

In order to understand the effect of the gravity of the 
central galaxy, ignored in the models described above, we 
have rerun model A3 including its contribution. As expected, 
being the galactic gravity dominant only in the inner ~ 20 
kpc, the basic results (cooling rate and variable profiles) 
are similar to run A3, and we do not show plots for this 
model. The only remarkable difference with model A3 is a 
reduced frequency of the jet, compensated with a larger av- 
erage power. At the end of the simulation the total energy 
injected by the outflows in these two models is almost iden- 
tical. 



3.14 Model A4, £= IQ-^ 

When the efficiency is increased further (e = 10~^, model 
A4, Fig. 4) the overall results are very similar to those for run 
A3. Again the jet heating generates fluctuations in the tem- 
perature which exceed the ones observed now in A 1795. The 
time averaged profiles (not shown), however, agree very well 
to the observations. At almost any time this model would 
be classified as a cool core cluster. 

The number of jet activations during the 7 Gyr of evo- 
lution is further reduced, with only a few happening within 
the first several Gyr. The total energy generated by the feed- 
back process is ~ 1.5 x lO"^ erg, again an acceptable value 
considering that the total 'available' BH energy is around 



3.1.6 Model A6. 



10" 



1.8 x 10*^^ erg {Ebk 



3.1.5 Model A5, 



O.lAfBHC^, with AfsH 



5 x 10" 



10^' 



This model (A5, e = 5 x 10~^, Fig. 5) clearly demonstrate 
the flaws of a too powerful feedback, albeit the total energy 
generated by the AGN is once more approximately the same 
as in all the other models. Now the 11 AGN outbursts oc- 
curring during the 7 Gyr evolution are quite violent, with 
power typically larger than 10*^ erg s~^ (right column of Fig. 
5). Notice as the outflow velocities approach the relativistic 
regime. As expected, the gas cooling is essentially zero, but 
the shock heating is too strong and the central temperature 
is too high for most of the time. With such an efficient feed- 
back cool core clusters would be a rarity, in contrast to the 
observational evidence. 



A feedback with e = 10"^ (an implausibly large efficiency, 
but such a model has been calculated for pedagogical rea- 
sons; model A6, not shown here) soon erases the initial cool 
core and the system would show a fiat temperature pro- 
file thereafter (apart an intermittent mini cool core, ~ 10 
kpc in size, associated with the central galaxy and expected 
also in non coo l core clust ers; iBri ghenti fc Mathewd |2002| : 
ISun et all [20071 : ISunll2009l '). In passing, we note that the 
ubiquitous presence of these galactic scale cool cores poses 
very strict constraints on the nature of the AGN feedback, 
which can not deposit a large amount of energy in the region 
near the central black hole. 



3.1.7 Summary of Feedback A (AAfcooij 

In this section we have illustrated the global features of feed- 
back A models, in order to check whether and when non 
relativistic outfiows are a tenable mechanism for AGN feed- 
back, able to shut off the gas cooling and at the same time 
preserving the 'cooling fiow' appearance {T{r) and n(r)). 
To summarize the main results, we find that the efficiency 
must be in the range 10~^ — 10~^ in order to generate suc- 
cessful models. Outflows must be relatively infrequent and 
powerful. The total energy necessary to (almost) halt gas 
cooling for 7 Gyr of evolution is about 1.5 x 10*'^ erg, not 
surprisingly of the same order to the energy radiated away. 
However, the energetic balance feedback energy « energy ra- 
diated away does not guarantee the success of a particular 
model (see models Al and A5, for instance). Instead, it is 
crucial the way this energy is communicated to the ICM, 
with the appropriate time-scales and power. 



3.2 Feedback B (intermittent) and C (continuous) 

We now illustrate the results when the jets are intermittently 
activated at predetermined times (feedback scheme B) or by 
assuming a continuous outflow (feedback C). 

In run Bl, outflows are generated every 200 Myr by set- 
ting iijct = 10^ cm s~^ inside the active region, for a time of 
20 Myr each. The results are shown in Fig. 6, rightmost col- 
umn. The general features of this model are similar to those 
of run A2. The cooling rate is very low for the first ~ 3 Gyr 
and then increases slowly up to ~ 80 Mq yr~^ at t — 7 Gyr. 
The variable profiles remind those of a pure cooling fiow 
calculation, especially at late times. This model is accept- 
able for the first few Gyr. However, the central gas density 
slowly grows with time until radiative cooling prevails over 
heating, causing the cooling rate to surpass the threshold 
of acceptability. The total energy injected by the outfiows 
is > 2 X 10^'' erg, a much larger value than in successful 
simulations adopting feedback scheme A, like A3. This huge 
amount of energy can not come from a single black hole of 
typical mass ~ 10^ M©. 

Some of the problems of model Bl can be cured in- 
creasing the frequency of the jet events as shown in run B2 
(second column from the right in Fig. 6). This model gener- 
ates outfiows every 100 Myr, each one lasting 10 Myr. Their 
velocity is the same as in run Bl. Although the gas cooling 
rate is much reduced up to i = 5 Gyr, it rapidly increases 
in the last ~ 2 Gyr. The temperature and density profile 
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Figure 4. Evolution of model A4 (e = 10 ^). Plots description analogous to Fig. 2. 
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gradients become too steep after ~ 1 — 2 Gyr, turning very 
similar to those of model CF. While these global properties 
would make this model acceptable, the required total energy 
is again ~ 2 x 10®^ erg, a value difficult to justify. 

A better model can be obtained increasing the fre- 
quency of the AGN feedback even further. Model B3, second 
column from the left in Fig. 6, illustrates a run with outflows 
activated every lO'^ yr (similar to the ou tbursts frequency es- 
timated for Perseus: iFabian et al.ll2006l '). each one persisting 
for 10® yr. As in models Bl and B2, the velocity is set to 
10'' km s~^. Only a very small amount of gas cools, and the 
variable profiles resemble again those of a classical cooling 
flow. The total feedback energy, ~ 1.6 x 10®"^ erg, is still 
inconveniently large. 

The left column of Fig. 6 shows the outcome of a sim- 
ulation in which the feedback is always active (scheme C). 
Here the outflow velocity is Ujot = 6000 km s~^ (model C2). 
This is clearly an extreme model that we have calculated 
mainly for pedagogical reasons. The ICM in the centre is 
continuously heated and transported to large radii, and al- 
most no gas is able to cool. The central density rises secularly 
and it is possible that some cooling would happen had we 
evolved the cluster for a longer time. The temperature pro- 
files agree very well with observations, with a flattening due 
to shock heating visible in the inner ~ 10 kpc. Another un- 
pleasant feature of model C2, besides the excessive injected 
energy, is that the continuous outflow forms a tunnel in the 
ICM where it propagates without inflating any cavity (see 
also the discussion of the Bondi accretion models in Section 
3.4). 

We have experimented the effect of varying the outflow 



velocity in feedback C simulations (not shown here). Increas- 
ing the jet velocity to Wjot = 10* km s~^ (model C3) leads to 
negative temperature gradient in the inner ~ 10 kpc, with 
r(0) ~ 5x10 ^ K, which may be in contrast with the observa- 
tions of lSunI (|2009). Conversely, with a reduced jet velocity 
iijct = 2 X 10^ km s~^ (model CI, also not shown here) the 
gas cooling can not be halted and this cluster resembles a 
pure cooling flow. 

In summary, it is possible to find partially successful 
models using feedback B or (especially) C, but it is difficult 
to justify the large amount of outflows energy required or 
the extreme character of a continuous jet, when we clearly 
see the signature of intermittent jets and AGN feedback in 
the different generation s of X-ray cavities in many clusters 
(McNamara fc Nulsen i2007l : I Wise et al.ll2007l : iFabian et all 
l2000h . 



3.3 Role of the jet size 

Finally, we address how the results depend on the size of the 
jet active region. We increased its length to ~ 17 kpc (the 
width being the same as before, 2.7 kpc) and calculated sev- 
eral models varying the efficiency. We do not show figures 
for these runs, but only briefly discuss their essential prop- 
erties. Overall, we find that the results are very similar to 
those described in the previous section. When the efficiency 
is e = 10~^ (A2L model), we find that the cooling rate grows 
steadily with time, reaching Mcooi ~ 80 M© yr~^ f = 7 Gyr. 
The temperature and density profiles resemble those for run 
A2, although they are slightly smoother for r <, 2Q kpc, as 
expected given that the outfiow shocks are generated at the 
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tip of the source region, located at z ~ 17 kpc. The run with 
e = 5 X 10""^ (A3L) has excellent attributes, with very low 
Mcooi ;$ 10 M0 yr"'^ at any time, and smooth profiles in very 
good agreement with those observed for A 1795. Finally, we 
find that high efficiency e = 10~^ still generates an excellent 
model (A4L), with very little cooling and very good density 
and temperature profiles, where a small temperature bump 
at r ~ 30 kpc, due to the jet shock, is sometime visible in 
the mass weighted, azimuthally averaged profile. 

At the other extreme, we run a model (ASS) in which 
the jets are generated by imposing a mass, momentum and 
kinetic energy flux in a small region Ax = Ay = 2.7 kpc 
at the centre of the boundary plane z = 0. Again the out- 
flows are triggered when gas cools, with a given efficiency. 
The velocity of the inflow is set by the requirement that 
the kinetic energy injected in a given timestep is eAMcooic^ 
(see Section 2.2). Th is sche me, si milar to tha t adop ted by 
IVernaleo fc Revnoldsl l|2006f ) and iHeinz eTall l|2006l ). does 
not directly change the ICM velocity, but the hot gas is 
pushed from below by the outflow, which enters the grid at 
z = 0. With efficiency e = 5 x 10"'^, the same as model 
A3, we find that the cooling rate (not shown) is acceptable 
(Afcooi ;$ 10 M© yr~^) until t ~ 5.5 Gyr, then increases up 
to Afcooi ~ 60 M0 yr~^ at the end of the calculation. The 
temperature profiles often show a positive central gradient 
in agreement with the observations of A 1795, although a 
sharp central temperature peak in the inner ~ 10 kpc is 
present for ~ 20% of the time. We did not pursue the search 
of the best model using this outflow generation method, by 
flue tuning the efficiency or other parameters. 

We conclude that the size of source region is not a key 
parameter of our models when the feedback power is linked 
to the cooling gas, as in the runs presented in Section 3.1. 



3.4 Feedback triggered by Bondi accretion 

As discussed previously, the feedback schemes adopted for 
the models presented in Section 3 are just a convenient way 
to link the response of the black hole to the cooling of the 
ICM, in order to make the feedback process self-regulating. 
The physics of the accretion, black hole growth and outflow 
generation, poorly understood in general, is by no means 
captured in these simulations. This is a general weakness of 
current models of AGN feedback: the energy (or momen- 
tum) is injected in the ICM according to essentially ad hoc 
prescriptions. 

A manifestation of this deficiency is the excessive ac- 
cretion rates occurring even in our most successful models 
using feedback A. For instance, in model A4 the most power- 
ful jets (Pjct ~ 10** erg s~^) imply accretion rates of ^ 2000 
Mq yr~^, much larger than the Eddington rate for a black 
hole mass of 3 x 10^ M©, MEdd ~ 60 M© yr~\ These s uper- 
Eddington values are not very common for AGN (e.g. iKina 
I2OO9I ) and yet, according to the results based on feedback 
A (Section 3.1), we need fast and energetic outfiows, which 
imply accretion of large gas masses, to shut down the gas 
cooling. 

At this point we tried a different approach to mitigate 
this too exp losive behaviour. The Bondi accretion theory 
(|Bondilll952l ). although highly idealized with respect to the 
complexity expected in real accretion systems (thin or thick 
discs, ADAF, etc.), has been widely used to estimate accre- 



tion rates on supermassive black holes. However, in presence 
of a standard cooling flow the 'unpertubed' gas is not at rest 
as in Bondi theory. The accretion rate on the central black 
hole is then determined by the ICM inflow rate and should 
be proportional to the cooling rate, which in turn is deter- 
mined by the global physical conditions of the ISM or ICM. 
Thus, the Bondi rate should not be particularly relevant if 
the gas cooling rate is large. 

Despite this theoretical consideration, our resolution is 
still very far from capturing the real central accretion en- 
gine (few tens of Schwarzschild radius) and we have to base 
the triggering mechanism to mean large scale values. In this 
sense it is interesting to calculate a model in which the AGN 
feedback is linked to the accretion rate estimated by the 
Bondi's formula: Mb = 47r(GMBH)^po/Cso (see Section 2.2). 
Note that the accretion increases when temperature drops 
and density grows, or better when the entropy (s oc T/p^' '') 
is falling off. Therefore, also Bondi models are sensitive to 
gas cooling, but in a more gentle and moderate manner. If 
such a small accretion rate is sufficient to halt the cooling, 
then this model is self-consistent. 

In a s a mple of nine X-ray bright elliptical galaxies 
lAUen et al.l (|2006l ) estimated typical Bondi accretion rates 
Mb ~ 0.01 Mq yr~^. For galaxies at the centre of rich clus- 
ters the Bondi rates are not expected to greatly exceed these 
values. Therefore, the Bondi AGN feedback model should 
operate in an opposite regime than that in models A3 or 
A4. In particular, we expect a quasi-continuous AGN activ- 
ity of moder ate intensity. Such a model is similar to that 
calculated bv lCattaneo fc Tevssied l|2007r ). who were able to 
drastically reduce the cooling rate (injecting also much ther- 
mal energy), although the price to pay for this result was 
the absence of a cool core. 

In the first BONDI model the outflow is generated in the 
same way as for feedback scheme A, that is injecting the en- 
ergy to the gas contained in the active region (2.7 x 2.7 x 6.8 
kpc''), but with an efficiency of 0.1. Note that here we aver- 
aged po and Cso in a radius of 10 kpc from the centre. We do 
not show figures for this run, but limit ourself to a general 
description. As anticipated, this model is characterized by an 
almost continuous activity of relatively weak outflows, with 
typical velocity of ~ 1000 km s"'^. The density and tempera- 
ture profiles are very similar to those of the classical cooling 
flow simulation and the cooling rate rises to unacceptable 
values of ~ 250 Mq yr~*. This simulation recalls Al: the 
very frequent, low power outflows are not able to prevent 
massive cooling, resulting in a model almost undistinguish- 
able from a pure cooling flow. It is important to note, that 
even if the heating is always turned on, the jet power is 
greatly oscillating from 10*^ to few 10*^ erg s~^. In fact the 
Bondi rate feels the low entropy gas which is cyclically build- 
ing in a small torus, around the perpendicular outflow (~ 10 
kpc). In the end this simulation fails to meet our requisites 
for a successful model. 

In order to partly detach the Bondi accretion rate from 
the huge amount of cooled gas in the inner tens of kpc (and 
thus Feedback type A), we need to double the resolution 
and then halve the averaging zone (5 kpc) . To explore Bondi 
feedback further, we wanted also to use the smallest jet pos- 
sible, calculating a second model with the alternative jet 
generation scheme: now the jet (mass, momentum, and me- 
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Figure 7. Evolution of model BONDI2. Here the jet is injected from the grid boundary at 2 = (see Section 2.2 for details). In the left 
column are represented the same quantities as in Fig. 2. In the right column are shown, from top to bottom, the total energy injected 
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chanical energy flux given by O.IMbc^) is injected from the 
grid boundary at z = as described in Section 2.2. 

The results for this model (BONDI2) are shown in Fig. 
7. The T and n profiles and the cooling rate are good, with 
the latter exceeding the observational limits only after ~ 6 
Gyr. As in the previous Bondi model the activity of the AGN 
is continuous and moderately weak, but now the jet power 
is very steady (on the contrary of previous simulation) in 
the range 2 — 3 x 10^* erg s~^, for a total injected energy 
of ~ 7 X 10®^ erg. In fact, the accretion rate does not vary 
much, Mb ~ 4-8 x 10~ ^ Mr y r~\ slightly larger than the 
estimates bv I Allen et all l|2006l ). 

In shrinking the central averaging zone (for the entropy) 
the Bondi accretion does not present many fluctuations and 
the triggered outflows seem capable of reducing the cooling 
flow. This might be more evident with higher resolution (up 
to Bondi radius) and thinner jets. 

Furthermore, comparing both Bondi models, we have 
noticed that, from a numerical point of view, injection 
through the boundaries seems a better method to produce 
a jet an d subsequent entrai nment (as pointed out by other 
authors, lOmma et al.l 12004 for example). This way the jet 
flux is tightly coupled to the (PPM) hydrodynamical algo- 
rithm and not inserted like a split source term, modifying 
by-hand flow variables in cells of the effective computational 
domain. Moreover, keeping a fixed jet density, the velocity 
can not reach very high values (like 10^ km s~^), facilitating 
the stability of the code and moderate CFL number (~ 0.4). 

A possible riddle for BONDI2 model is the absence of 



frequent jet-inflated spherical cavities (see also model C in 
Section 3.2). The continuous AGN activity (mean v ~ 7000 
km s^^) carves a narrow tunnel of about 50 kpc in length, 
although its density contrast with the environment is large 
only for z < 20 kpc (and not particularly evident in the SB 
maps). See Section 6.3 for a deep discussion. 

We will certainly extend the study of Bondi-type feed- 
back in a dedicated future paper. 



4 DYNAMICS OF MODEL A3 

In this Section we will discuss the dynamical evolution of the 
outflows for one of the best models, A3 (e = 5 x 10"'^). Con- 
trary to the almost steady, and not variegate, evolution of 
model BONDI2, the density and temperature maps (i.e. the 
X — z midplane) of A3 show signiflcant temporal variations. 
Note that in this Section we are analysing physical quanti- 
ties, while the right comparison with observations should be 
done through emission-weighted ones. We will broaden this 
aspect in Section 5. 

The flrst outflows is triggered at t ~ 160 Myr (see Fig. 
[3|, with a velocity Hjet ~ 2 x 10* km s~^, mechanical power 
Pjct ~ 1.3 X 10** erg s~^ and a total energy injected of 
~ 4.5 X 10^° erg. The large jet power is due to our adopted 
feedback method A, where outflows last only while gas is 
cooling to low temperatures. This implies that the AGN 
is active only for few timesteps (typically At ~ 10^ — l(f 
yr), since the jet promptly stops the gas cooling in the cen- 
tral region. During a single timestep large masses of gas can 
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Figure 10. A3 model: mass weighted, azimuthally averaged tem- 
perature profiles for the six times showed in upper two rows of 
Fig. 9. 



cool (also a result of the coarse resolution used), therefore 
the energy eAMcooic^, corresponding to the cooled gas mass 
AMcooi, results in a large, and probably exaggerated, in- 
stantaneous outflow power. 

In Fig. 8 and 9 are represented density and temperature 
maps in the x — z plane (the outflows are generated along 
the z— axis) for the first two AGN outbursts. At i = 170 
Myr (Fig. 8a), after ~ 10 Myr since the jet started, an el- 
lipsoidal cavity with major (minor) semiaxis of about 15 
(10) kpc has been carved, surrounded by a shell of shocked 
gas, whose density is about twice the ambient density at 
nearby locations. The temperature of the shocked gas is 
~ 10* K, about three times that of the unperturbed ICM. 
Therefore the young cavity, expanding approximately at the 
sound speed is surrounded by a weak, hot rim. The cavity 
has a relatively low|j density contrast with the environment, 
~ 3 — 5. Outside the cavity shock the ICM is still slowly flow- 
ing in, as a classical cooling flow. The azimuthally averaged 
(mass weighted) temperature profile (Fig. 10) shows a spike 
at r ~ 30 kpc, with a fractional increase ~ 60%. The cavity 
expands, increasing its ellipticity, and at i ~ 180 Myr (Fig. 
8b) it has approximately a cylindrical shape, extending up 
to ~ 100 kpc and with a radius --^ 15 kpc. The bubble is 
filled with hot (> 2 x 10* K) gas raising along the z-axis at 
about its sound speed. The density contrast is again 3 — 4. 
A very weak 'pear-shaped' shock surrounds the cavity. The 
density gradient along the post-shock region, with denser gas 
closer to the equatorial x — y plane, makes the low— z part of 
the shock detectable, while at large distance from the centre 
the low density contrast likely prevents the shock detection 
(see also the X-ray brightness map in Fig. 11). The physical 
temperature jump across the shock for 2 < 25 kpc is only 
~ 25% and increases slightly at large z. The structure qual- 
itatively re minds the one observed in the elliptical ga laxy 
NGC 4636 l|Finoguenov fc JonesHJOoH : iBaldi et al.ll2009l ). In 
the azimuthally averaged temperature (Fig. 10) the weak 
shock and the heated gas are visible as a small bump ~ 25% 
in amplitude, located at 10 < r < 70 kpc. 

^ We lack the effect of relativistic protons, whose pres- 
sure is likely to be rel ev ant in e xpanding the cavity 
llMathews fc Briehentill2008allbt IGuo fc Mathewsii20ia) . 



20 Myr later the cavity lengthens and narrows, while the 
back How generates a dense, relatively cold filament protrud- 
ing in the cavity (e.g. [Mathews & Brighenti 2008a; Gar dinil 
[20071) . As for the previous times, in most of the volume of the 
cluster the ICM is inflowing as in a standard cooling flow, 
and the averaged profiles agree very well with those observed 
for A 1795; after 40 Myr since the powerful AGN outburst, 
the cluster fully restored the 'cool core' appearance it had 
at the beginning of the calculation. A very small amplitude 
'ripple' {AT/T ~ 10%) is visible in the temperature profile 
at ~ 100 kpc as the integrated effect of the elongated weak 
shock. 

At t = 250 Myr (Fig. 8c) the filament reaches 2 = 100 
kpc, while the cylindrical (subsonic) outfiows in the (now 
almost disappeared) cavity reached z ~ 250 kpc. The tem- 
perature ripple moves forward at the sound speed and slowly 
weakens. 

At i ~ 350 Myr (Fig. 8d) the cavity disappears while 
the dense, relatively cold filament formed by gas formerly 
at the centre and lifted at large 2 by the jet motion, is still 
clearly visible. The inner part of the filament, at 2 < 40 kpc, 
reverses its velocity and starts to fall back toward the centre. 
This is reminiscent of the kinematic o f the emission lin e 
filaments observed in the Perseus cluster (IHatch et al.ll2006l ) . 
The outer region around the 2— axis is instead still slowly 
fiowing out. This moment marks a new phase in the cluster 
lifecycle. The fall back of relatively dense gas preludes the 
next cooling/feedback event, which occurs at t ~ 450 Myr. 
The feedback cycle starts again in a qualitatively similar 
way. 

Fig. 8e shows the density map at 465 Myr. A new cavity 
is formed, at first of approximately spherical shape centred 
at 2 ~ 25 kpc, radius ~ 20 kpc and high density contrast 
(50 — 80). The almost spherical symmetry of the cavity is 
caused by the infiow along the jet axis, whose ram pressure 
slows down the expansion along that direction. The shock, 
already weak (Mach number ~ 1.5) is also nearly spherical, 
with radius '-^ 35 kpc. Note that, at this time, the outfiow 
has an energy of ^ 8 x 10®" erg, while the cavity has roughly 
~ 5 X lO'^^ erg (with the usual 2.5 PV). Thus most of the 
mechanical energy is not used to form the cavity. 

Again the temperature profile shows the signature of 
the young AGN outburst with a strong peak for r < 50 kpc, 
approximately the location of the shock along the 2— axis. 
A very weak ripple, vestige of the first jet event, is visible 
at r ~ 400 kpc. 

At i = 500 Myr (Fig. 8f ) the cavity shape is strongly af- 
fected by the back flow and the density contrast lowers. The 
weak shock (M ~ 1), slightly elongated in shape, is now 
located at a distance of ^ 100 kpc, much further away than 
the cavity. Within r ~ 100 kpc the cluster atmosphere is 
slowly moving outward (with velocity in the range 200 — 500 
km s~^). Being the motion very subsonic the density pro- 
file changes very slowly. The outfiow is decelerating and in 
~ 100 Myr reverses its direction approaching the dynamics 
of a classical cooling fiow. Thus, the ICM in the cluster core 
undergoes cycles of slow contraction and expansion, follow- 
ing the rhythm of the AGN activity. 

The cylindrical outfiow with average velocity Uz ~ 300 
km s~^ along the 2— axis, remnant of the previous AGN 
outburst is still present, and extends beyond 2 — 200 kpc. 
At the same time the infiowing gas in the lower part of the 
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Figure 8. A3 model: logarithm maps of electron number density (cm '^) in the x—z midplane (kpc unit), with velocity field superimposed. 
The color scale is given by the bar on the top, while arrows length normalization varies. Times are indicated on every top right corner. 
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Figure 9. A3 model: maps of the logarithm of gas temperature (K). See Fig. 8 for other details. 
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Figure 11. X-ray surface brightness maps for model A3 at various times. The a:;-axis is horizontal, while the 2;-axis is vertical (kpc units). 
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Figure 12. Emission weighted temperature maps for model A3 at various times (see also Fig. 11). 
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filament is effectively preventing the cavity from expanding 
or raising buoyantly above 2; ~ 50 kpc. 

Finally, in the bottom row of Fig. 8 we show few snap- 
shots of the density distribution at late times. As the cycle of 
feedback proceeds the flow develops a more turbulent char- 
acter. In panel 8g is represented the density map at i = 3 
Gyr. The cavity, generated at t ~ 2.91 Gyr, is distorted by 
the ascending backflow and by falling gas in the filament 
along the z— axis. As a result it acquires the shape of an 
asymmetrical torus. As in the previous aftermaths of the jet 
episodes, the ICM is fiowing inward in most of the cluster 
volume, and is deviated outward along the z— direction when 
it reaches r ~ 20 kpc. 

Panel 8h illustrates the cluster during a quiescent period 
(t = 3.5 Gyr), just before a jet is triggered. The density 
distribution is very smooth and the cluster appearance is 
that of a standard cooling flow (see also the profiles in Fig. 
3). A large, slightly overdense region around the z— axis is 
falling toward the centre with velocity varying between 50 
and 500 km s~^. 

Finally, in panel 8i we show the density distribution at 
t = 6.5 Gyr, again in an epoch long after an AGN outburst. 
Again, the density is smooth in the core region, while shows 
some variation for r > 50 kpc. On the contrary, the velocity 
field is chaotic and very subsonic; it promotes mixing of the 
metals produced by the SNIa exploding in the central galaxy 
(see Section 5.3). Few streams of moderately overdense gas 
are falling from large radii r ^ 100 kpc with velocity ~ 300 
km s^^. 
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Figure 13. A3 model: X-ray surface brightness and emission 
weighted temperature ID cuts (of Fig. 11 and 12 maps), through 
2 = 25, 50, 50 kpc, at 465, 500, 3000 Myr, respectively. 



5 X-RAY OBSERVABLE PREDICTIONS 

The x — z cuts through the centre of the cluster are very use- 
ful tools to investigate the intrinsic dynamical evolution of 
fiow variables, like density, temperature and velocity. How- 
ever, in astrophysics we are limited to observations based on 
the surface brightness (SB) in a typical spectral range (in our 
case mainly X-ray band). In theory we have to mock every 
aspect of the observation, from different atomic emissions 
to instrument response. Due to our limited resolution, it is 
sufficient for us to just integrate the emissivity or emission- 
weighted quantity along line of sight (j/-direction) in an en- 
ergy range of the X-ray band ~ 0.5 — 10 keV (similar to 
Chandra) . Note that we have created a parallel code that is 
able to interpolate every slice of the data cube (remember 
that we have 9 different levels) and then perform the above 
mentioned integration, in order to obtain more precise SB 
maps. With these maps we can test further the observability 
of our models and, in the lack of real ones, just make pre- 
dictions that could be in the near future verified or falsified. 
We will also test common assumptions regarding gas hy- 
drostatic equilibrium, taken in many observational analysis, 
that seem too restrictive. 



5.1 Cavities and shocks 

We begin investigating the detectability of faint features, like 
X-ray cavities and shock waves, generated by the AGN out- 
flows. In Fig. 11 we show the X-ray surface brightness maps 
for run A3, at the same times showed in Fig. 8 and 9. The 
emission weighted temperature map is displayed in Fig. 12. 



X-ray depressions are clearly seen at t = 180, 465, 500 and 
3000 Myr. Evidently, subrelativistic, massive, coUimated 
outflows can generate cavities with typical diameters of 15- 
40 kpc. Often the X-ray bubbles are surrounded by bright 
rims of shocked gas. Relatively weak waves are present at 
large distance from old cavities, both features being gener- 
ated by the same AGN outburst. 

To better quantify the brightness perturbations pro- 
duced by the outflows we show in Fig. 13 the profiles along 
the a;— direction at t = 465, 500 and 3000 Myr, taken 
at z = 25, 50, 50 kpc, respectively. The sharp jumps at 
X ~ ±30 kpc at t = 465 Myr, where the cavity shock in- 
creases the brightness by a factor ~ 2.5 would be manifestly 
visible in the X-ray image. The relative central depression 
(the cavity) is indeed brighter than the same region before 
the jet activity. The emission weighted temperature map re- 
veals that the cavity region is markedly hotter (by ~ 80 
%) than the nearby ICM, which is not commonly observed. 
This is a consequence of the infrequent, but explosive out- 
flows like those of model A3. 

The maps and profiles at t = 500 Myr show a weak cav- 
ity centred in y '^ 30 kpc, with radius ~ 20 kpc. This cavity 
is also slightly hotter than the surrounding gas. At this time 
the most interesting feature is the weak shock located at a 
radius ~ 100 kpc: Mach number is ~ 1.1, in good agreement 
with observations of shocks driven by AGN ( Blanton et al.l 
120091 ). The surface brightness profile shows a clear front at 
the shock position and the emission weighted temperature 
jumps of ~ 10 %. 

At i = 3 Gyr another weak cavity, again ~ 20 kpc in 
radius, is present. The brightness depression is only ~ 10 % 
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Figure 14. Emission-weighted iron abundance maps {Zq unit) in the x ~ z midplane for model A3, at three different late times. 



and the temperature is higher than the ambient one by the 
same percentage. 

To summarize, our best models presented in Section 
3.1, with relatively powerful and infrequent outflows, have 
the tendency to produce cavities in a violent way, generating 
shocks which heat much the surrounding ICM. This is also 
a reason why we investigated the effect of a weaker and 
quasi-continuous self-regulated feedback (Section 3.4). We 
will discuss the differences in Section 6. 



5.2 Iron enrichment and mixing 

It is expected that directional outflows would generate metal 
inhomogeneities through the ICM. Iron is the most rele- 
vant and easily measured element, being produced mostly 
by SNIa, which are still exploding in the giant elliptical at 
the centre of cool core clusters. The same outflows, however, 
generate turbulence and bulk motion, which in turn tend to 
stir and mix the ICM, restoring homogeneity and erasing 
abundance gradients. 

We present here a brief analysis of the emission- 
weighted iron abundance evolution for model A3. We model 
only the Fe-enrichment produced by the SNIa (and stellar 
winds) occurring in the central galaxy, in the time interval 
6.7 — 13.7 Gyr, the latter being the current age of the uni- 
verse tn. Neglecting the iron produced by the SNII and by 
the other cluster galaxies, we are not in the position to in- 
vestigate the complete chemical evolution of the cluster (this 
will be the subject of a future work). Instead, we are inter- 
ested here in quantifying the abundance anisotropies caused 
by the AGN outbursts. 

The current SNIa rate is assumed to be ~ 0.1 SNu (su- 
pernovae in 100 yr per lO^^Ls,©), sl ightly below to that es - 
timated in local early-type galaxies (|Cappellaro et al.lll999l ) 
but in agreement with the value necessary to generate the 
obs erved abundance in well observed giant elliptical galax- 
ies (|Humphrev fc BuotdlioO^ : iMathews fc Brighentill2003h . 
The time de pendence of t he SNIa rate is assumed to be 
« (i/i„)-'-' llGreggid l2005l'l. 

The iron density is implemented in the code as a tracer 
of the flow (ppc — 4>Ycp with ^pe a scalar between and 1), 
following the usual advection equation: 



dpt 



dt 



+ V • (pfcv) = Spo 



where the source term Sfc dep ends mainly on asNp* (Sec. 
2), as previously discussed (see lMathews fc Brighentill2003l 
for more details). 



In Fig. 14 we present the Zfc maps for model A3, at 
three different late times. Note that the background (black) 
is zero, but in reality should be around 0.3 Zq: here we are 
just interested in the contrast. At 3 Gyr (first panel) the 
abundance is highly asymmetric. A powerful jet event has 
recently occurred (see Fig. 3) , therefore the metals created in 
the central elliptical galaxy (within few kpc) are transported 
out along the 2;— direction up to a distance of 150 — 200 kpc. 
This is an unmistakable mark of the presence of a powerful 
AGN outflow. In fact we suggest that iron abundance can 
be used as a reliable tracer for any AGN jet-outflow activity, 
instead of common entropy maps. At the centre the abun- 
dance is about 0.4-0.6 Zq, while further away along jet axis 
the contrast is 0.1-0.2. This kind of featu re is supported by 
recent deep Chandra observations done bv lKirkpatrick et al.l 
(2009). They produced Fe-maps showing the enrichment of 
gas along the jets of Hydra A, a quite massive galaxy cluster. 
It is striking that our simulated maps resemble these obser- 
vations also in a quantitative way. Other authors, Doria et 
al. (in preparation), found a similar behaviour in a different 
cluster dominated by AGN activity (RBS797). 

At 3.5 Gyr (second panel) we are in a period of rela- 
tively quiescence and so the now dominating (AGN induced) 
turbulence and vorticity promote the iron diffusion. At the 
centre, the new iron cloud associated to the cD galaxy is 
clearly detached from the older ejected material, the latter 
covering now a more uniform area. The iron diffusion is more 
evident at 6.5 Gyr (third panel), again in a moment of AGN 
quiet, where the iron is very diffuse, with a low enhancement 
of 0.05-0.1 Zq. 

In the end we can affirm that in a single cycle the iron 
abundance passes from a phase of high asymmetry along jet 
axis, when the outflow has recently turned on, to a phase 
of turbulence and mixing, in which spherical symmetry is 
almost restored. The iron gets spread within few 100 kpc, 
just as observed: without the influence of an AGN such ma- 
terial would be indeed confined near the effective radius of 
the central galaxy. 

5.3 Hydrostatic equilibrium 

In this Section we investigate the effect of the perturba- 
tions generated by the AGN outflows on the mass deter- 
mination using X-ray observations (through Eq. (O). A 
standard method to estimate galaxy cluster mass profiles 
is to assume spherical symmetry and estimate the gravita- 
tional mass using the hydrostatic e quilibrium equation (e.g. 
iBuote et al.ll2007l : IVikhlinin et al.li2006 for two recent com- 
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pilations). When the flow velocity is much lower than the 
sound speed, the hydrostatic assumption is fully justified. 

It is well known that in real clusters turbulence or large 
scale m otion can systematically bias the mass estimate by 
<, 20% (iRasia et alJlJOO ^: 'Nagai et al.ll2007l : iMahdavi etlo] 
l2008l : |Pifraretti fc Valda rnini 200i). Therefore, it is interest- 
ing to investigate how much the flow perturbations gener- 
ated by the AGN feedback affect the mass measure. 

Here we quantify the discrepancy between the estimated 
gravitational mass and the real one for a model adopting 
feedback scheme A and one assuming the more quiet feed- 
back Bondi. We have calculated the azimuthally averaged 
density and emission weighted temperature profiles in rings 
of progressively increasing width, from 10 kpc in the cen- 
tre up to 100 kpc at large radii. These profiles, inserted 
in the hydrostatic equilibrium equation give the estimated 
mass profile. This procedure is not perfectly accurate: we 
should have 'observ ed' our simulations with specific soft - 
wares li ke X-MAS (iGardini et al.l |2004| : iRasia et al.l 120081 ) 
or XIM IjHeinz fc BrueggenI 12003 ) to properly compute the 
averaged profiles. However, it is not the purpose of this pa- 
per to thoroughly investigate this topic and our approximate 
analysis is sufficient to make our point. 

The calculated mass profiles at various times, as well as 
the exact mass profile, are shown in Fig. 15 for model A3 
(top panel) and BONDI2 (bottom panel). 

In run A3 the powerful outflows, although not able to 
significantly perturb the thermal state of the cluster (which 
always preserves the cool core appearance), are effective in 
disturbing the quasi-hydrostatic equilibrium present in the 
pure cooling flow model. This results in a typical error in the 
mass determination of a factor 2-3. Most likely the estimated 
mass is lower than the real one, and the discrepancy is larger 
in the region r < 100 kpc, where the feedback affect the ICM 
the most. 

In model B0NDI2, instead, the outflows are 3-4 orders 
of magnitude less powerful and, despite the fact that they are 
continuously generated, they seem rather innocuous for the 
dynamical state of the ICM and the hydrostatic equilibrium 
approximation is safe. The computed mass proflle agrees 
very well with the real one, with slight discrepancy only 
visible for r < 20 kpc, a region inadequately resolved in our 
simulations. Of course, the weak, continuous jets present 
in this run are incapable to generate cavities (see Section 
3.5), hence it is likely that real clusters undergo at least few 
stronger AGN outbursts, after which the dynamics of the 
ICM would be similar to that described for model A3. 

In summary, we expect that for clusters where X-ray 
cavities and/or shocks are present the total mass estimated 
through the assumption of hydrostatic equilibrium might 
be in error by a factor of ~ 1.5, occasionally the error could 
get as large as a factor of 2-3. Conversely, with gentler con- 
tinuous outflows generated by the relatively accretion rate 
predicted by Bondi theory, the estimated and the real mass 
are in excellent agreement, with errors always below 8% (see 
also iGuo fc Mathews! |2010| ) . 



6 DISCUSSION 

In this paper we have presented several moderate resolution 
simulations of the interaction between AGN outflows and 
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Figure 15. Mass profiles for model A3 (top panel) and BONDI2 
(bottom panel), caleulated using the hydrostatic equilibrium 
equation at times separated by 0.5 Gyr (model A3) and by 1 
Gyr (model BONDI2), and the exact mass profile (red dashed 
line). 



the ICM in massive galaxy clusters. The purpose of these 
calculations was to investigate if a purely mechanical AGN 
feedback of this kind is able to solve the so called cooling 
flow problem: the dearth of cooling gas in cool core clusters 
with short central cooling times. 

The necessity of covering a large range of spatial (kpc 
up to 5 Mpc) and temporal scales (fraction of Myr up to 7 
Gyr) limited our resolution of the very inner accretion re- 
gion. Hence we had to link the feedback to some large scale 
mean quantity, like AMcom or MBondi, with the obvious re- 
sult of some discrepancies with observables (like bubbles) . In 
any case gas accretion onto the SMBH and subsequent jet- 
outflow ignition is still obscure, in particular the amount of 
ICM entrained and shocked. Furthermore, we do not have 
a long term (Gyr) evolution of theoretical models for ac- 
cretion, due to computational limitations. In the end our 
simplifled scheme seems a good and efficient method to test 
AGN feedback on galaxy cluster scales. 



6.1 Comparison with other works 



In the field of simulated 
the w ork done in the 
20021 : iBrighenti fc M athewj 
2004: IPall a Vccchia et al 
Siiacki et al.l 



AGN feedback most of 
past (IBrueggen fc Kaiser 



l2007l: 



J2003' 



2004: 



Ruszkowski et al 



Bruoggcn ct al. 2005 

Brueggen fc Scannapiecd 12009 



uegg _ _ ■ 

Mathews fc Brighentil l2008al lbl: iMathewsl |2009| ) has fo- 
cused on the creation of 'artificial' bubbles off-centre. This 
scenario is motivated by the fact that relativistic radio jets 
(distinct from our outflows) do not entrain much gas and 
just thermalize the ICM at the hotspot, generating a cavity. 
The difference with our study is that we employ a 
momentum-driven heating, instead of a pressure-driven one. 
As we have seen, purely kinetic feedback changes the whole 
dynamics. Intermittent bubbles are naturally created, but 
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this shocked and expanding gas is just one element of the 
heating process : the bow shock, the cocoon with entrained 
gas, the mixing through turbulence and vorticity are all fun- 
damental gears of the machine, everyone dominant at differ- 
ent phases of evolution. Moreover, artificial (hydrodynami- 
cal) bu bbles becomes unstable rapidly; on the contrary some 
studies l|Sternberg et al.ll2007l ) indicate that jet-inflated ones 
may stay intact for over 40 Myr, thanks to the vortex formed 
inside. It is probable that a combination of the two mech- 
anisms is required (outflow plus artificial buoyant bubble), 
even if it still unclear which of the two dominates. 

It is difficult to compare o ur work with other ki- 
netic outflows simulations (e.g. iRuszkowski et al. l2004l : 



lOmma et aLll2004l : IZanni et al.ll2005l : ISternberg et alll2007l) , 
because most of them do not implement radiative cooling 
and an initial cool core state. Additionally, they last for few 
hundreds Myr and are not suited to study the 'cooling flow 
problem', without the long term evolution. 

The only investigations done in this d irection are those 
of BM06 and ICattaneo fc Tevssierl l2007l . In the first one, 
as many previous studies, the jet velocity and power were 
set by-hand following observational estimates. We have seen 
(feedback scheme B), that in imposing such predetermined 
conditions, the mechanical feedback is not linked to the 
cooled mass or lower gas entropy (in quantity and time). 
The results are models with temperature and density pro- 
files si milar to a pure cooling flow o r with high cooling rates. 

In ICattaneo fc Tevssierl (|2007]) the injection scheme is 
self-regulated by Bondi accretion (for 12 Gyr), but after a 
weak burst of ~ 10*^ erg s~^ the profiles do not present the 
cool core appearance (flat temperature). Also our successful 
Bondi model is able to prevent the cooling catastrophe. The 
difference is that our profiles present ripples and signature of 
weak shocks. It is again difficult to make a full comparison, 
because a large fraction of their injected energy is thermal, 
not only kinetic. 

iDubois et al.l (|2010l ) expanded the previous model in a 
cosmological resimulation. Even if they focused on the BH 
growth history, they found Bondi outflows can prevent very 
peaked density profiles, despite the presence of negative tem- 
perature gradient at various times. Here the main difference 
with our simulations (and thus some results) is probably 
the cosmological context. We calculate a fully relaxed clus- 
ter, while they analyse a 1:1 merging event. In any case, we 
underline the fact that all these results give one clear indi- 
cation: AGN outflows are a key component of the feedback 
in galaxy clusters. 

The last comparis on is the one with 
IVernaleo fc Revnoldsl (|2006r ). They implemented a self- 
regulated mechanism similar to our boundary injection 
scheme. Their simulations present very different results 
from our (and Cattaneo et al.) analysis. After few hun- 
dreds Myr (with every AGN model) the cooling becomes 
catastrophic: > 2000 Mq yr~^. The problem of the creation 
of an unidirectional channel and subsequent deposition of 
energy is a difficulty found also in our outflow models, but 
partially resolved in the long term (see Section 6.3.1). 

In the next Sections we will deeply discuss the merits 



^ usually most of the mechanical energy is not used to form the 
cavity (see Sec. 4) 



and flaws, previously introduced, of our simulated feedback 
models. 



6.2 Cold-regulated explosive feedback (A) 

In all our study we covered a wide 'zoology' of triggering 
mechanisms, varying the mechanical efficiency and injection 
method. 

In the first series of models we linked the power of 
the AGN outflows to the cooling rate (feedback scheme 
A, Sections 3.1), the latter being calculated as AMcooi/Ai, 
where AMcooi is the mass of gas cooled within the region 
r < 10 kpc in a single timestep At. It is clear that the 
energy of a single jet event depends on numerical resolution 
through the Courant-Friedrichs-Lewy (CFL) stability condi- 
tion. Lower resolution calculations, with larger At and thus 
larger AMcooi, will generate more powerful outburst events. 

Although in the models described in Section 3.1 the 
total kinetic energy injected at a given time rests on the 
integrated cooled gas mass (which is only weakly depen- 
dent on the numerical resolution) , and can in principle vary 
among different models (i.e. different e), we find, somewhat 
surprisingly, that it is insensitive on the efficiency e and al- 
ways about 2 — 3 X 10^^ erg, a value comparable to the total 
energy radiated away. Evidently, the self-regulation mecha- 
nism assumed is very effective. 

The first astrophysically important result of these cal- 
culations is that the effect of the feedback crucially depends 
on the modality of the energy injection, not only on the total 
amount of energy. We have found that when the efficiency is 
low and the outflows are frequent and relatively weak (mod- 
els Al and A2), the gas cooling proceeds at a rate much 
higher than the limit allowed by observations. 

Conversely, when e ^ 0.01 (models A5 and A6) the 
cooling rate is reduced to a value fully consistent with the 
observational constraints. However the ICM in the core is 
heated too much and the temperature and density profiles 
do not resemble those of a typical cool core cluster. 

Only models A with efficiency in the range 5 x 10"'' ^ 
e < 0.01 are able to reduce the cooling rate to acceptable 
values preserving at the same time a central positive tem- 
perature gradient. Moreover, these models generate cavities 
and shocks, another important requirement that a plausible 
heating scenario must satisfy. Notice that also observational 
data (Merloni & Heinz 2008; La Franca ct al. 2010) suggest 
an average mechanical efficiency of AGN around 5 x 10""'. 

As pointed out before, other more ad hoc types of feed- 
back (not directly linked to cooled gas), like intermittent (B) 
or continuous with fixed velocity (C), present a cooling flow 
appearance, in the former case, or a too high total energy 
injection, in the latter. These models are not totally disas- 
trous and show some good features (as found in BM06) , but 
in the long term they are not satisfactory. 



6.2.1 Riddles 

Adopting a severe critical point of view, successful A mod- 
els may present two difficulties. First, gas often cools and 
accretes onto the black hole at a rate far above the Edding- 
ton li mit (sectio n 3.4), a fact which make them unpalat- 
able (JKinelbOOgf ). but not impossible (especially at higher 
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redshift, see for example ISzuszkiewca 12004 : iGhosh et alj 
l2010f ). It is known that BHs grow mainly through the 'QSO 
phase', that is throu gh radiative accretion at early times 
l|Hopkins et al.1 12006| ) . Therefore our accretion rates seem 
overestimated, considering that the BH mass should not in- 
crease much in the entire simulation time. 

Remember, however, that our resolution does not per- 
mit to track the exact amount of accreted material on sub-pc 
scale. Hence it is probable that only a small part of AMcooi 
falls onto black hole, and the other shall be just entrained 
or ejected by the outflow. It may be possible that the real 
efficiency could reach higher values (instead of 10~^ — 10"^ 
adopted in our successful simulations). 

For example, follow this consideration and let us define 
the inflow rate as the sum of the real BH accretion rate and 
the gas outflow rate: Min = Mace + Mout. If r; = Mout/Macc, 
then the jet power Pj<,t = treaiMaccC^ = £reaiMnC^/(l -I- ri). 
Thus we can say that the mechanical efficiency adopted in 
our simulations could be e = ercai/(l + rf). K mean r] can be 
retrieved as AMin/AMacc — 1- Now, assuming a small BH 
mass variation in 7 Gyr of 0.1 Mbh (that is AMacc — 3 x 10^ 
M0), and noting that for model A3 AMin = AMcoid ■^ 
3 X 10^ M0, then clearly 77 ~ 9. With this in mind, the real 
efflciency may be trcai = (1 + rfje ~ lOe = 5 x 10~^. Not 
only, we can also estimate AMout — 2.7 x 10^ M© and, with 
a duty cycle of 10% Q the mean outflow rate results in a few 
Mq yr~^. From here we can check the entrainment (or 'mass 
loading'): Aen = Mact/Mout, where Mact is the outflow rate 
in the active region. Taking a mean value of 10'^ M© yr~^, 
the entrainment is Aon ~ 10^ — 10'^. 

Regarding entrainment, in this flrst series of 3D simu- 
lations we just focused on the global consequences of mas- 
sive outflows on the thermal and dynamical ICM evolution 
on large scales. We assume that the outflow is momentum- 
driven with negligible thermal energy. The detailed process 
of generating the entrainment is here not a ddressed. One 
possible explanation is given by ISokeiJ |2008| . whose model 
forms slow massive wide jets. 

Summarizing, in order to have a small BH growth, in 
one of our best models, we have to assume massive out- 
flows with a rate ~ 9 times the accret ion one. The rati o r) is 
not well constrained, but for example iMoe et al.l (|2009l ) give 
values around 10. Nowadays nobody knows how much gas 
is really loa ded by the original jet, bu t ~ 10^ — 10'^ seems 
reasonable (|Cattaneo fc Tevssieii 120071 adopted 100 for ex- 
ample). This could be in the future another observational 
constraint for our numerical models. 

We conclude from all this reasoning that the high 'ac- 
cretion' rates in our successful simulations can be easily re- 
duced, introducing an higher (real) efficiency and consider- 
ing that only a fraction of AMcoid falls onto BH, while the 
other part is ejected and entrained. 

The second fiaw in models A may be that the process of 
cavity formation heats the surrounding ICM too much (Sec- 
tion 5.1): soon after jet ignition the temperature inside the 
bubble rises to high value (~ 10* K), and the ICM around 
is highly shocked. On the other side many observations 
llFabian et al.1 l200d: fMcNamara et al ." 2000; Blant on et al 



I2OO3I , I2OO6I : MN07) suggest that the bubble inflation should 
be gentler, in order to produce rims of gas cooler than the 
ambient medium. Again our resolution does not permit an 
accurate and detailed study of this and other kind of features 
(e.g. cold filaments). Certainly, a direct cause of the violent 
behaviour is that the injection energy in a single timestep 
is too much. This will probably be avoided with higher spa- 
tial resolution, because of smaller At and therefore smaller 
AMcoid and blowing instant power. Unfortunately with cur- 
rent computational resources we were limited to a very short 
evolution. In any case we plan to develop such a study in a 
future work. 

We can also argue that observations may be seeing bub- 
bles at a late time, well after their generation, and this could 
explain lower temperature and SB jumps. In fact of the three 
cavities analysed in Fig. 13, just the newly born one present 
sharp contrasts, while the other are older and hence feebler. 
Moreover, we can not exclude that in the entire simulation 
some bubbles are created in a gentler way by weaker bursts 
(~ lO''^ erg s-i). 



I2OOII : iNulsen et al.1 120021 : lBlanton"eral.. .2003 : .Fabian et al 



Thus a single jet event has AMout ~ 10^ - 10* M© 



6.3 Hot-regulated gentle feedback (Bondi) 

To circumvent some of the aforementioned difficulties en- 
countered with feedback A method, we changed completely 
direction, calculating models where the accretion rate was 
set to the Bondi rate. Contrary to models A, this generates a 
continuous feedback of moderate power (typically few x 10'*'' 
erg s~'). In order to be efficient, such non explosive jet power 
must be sufficiently steady in order to reduce cooling flow. 
As pointed out in Section 2.2, our resolution is far above 
Bondi radius {^ 50 pc), thus we have to rely again upon 
mean large scale quantities. In this case we must stay as 
much possible close to the BH, avoiding the simplest case 
of a few cells, because of numerical fluctuations. When the 
radius of (p, Cs) averaging is fairly small, ~ 5 kpc, the power 
is steady and the model successful: cooling rate are low, pro- 
files follow observations, and energies are contained. 

On the other side, with rav ~ 10 kpc, we revert again 
to a spasmodic feedback very sensitive to central cold mass. 
Hence this Bondi model recalls low efficiency Al model. This 
behaviour is also emphasized by the fact that here we are 
injecting mechanical energy directly in the domain cells, and 
not through a boundary (like BONDI2). In a short bursting 
event like models A this difference is not relevant, while with 
a continuous outflow, pushing the gas from below and not 
changing suddenly internal flow variables, seem to be a more 
efficient way to heat the inflowing central medium. 

This illustrates the (unfortunate) sensitivity of the sim- 
ulations on some numerical detail. The cause is that any kind 
of AGN feedback simulation with large spatial and tempo- 
ral scale integration must require a few quantities set by- 
hand. In general a fully self-consistent simulation is rather 
Utopian, for now. In our case the ignorance of detailed AGN 
accretion and jet-outflow physics, plus today computational 
resources, limits us to simplified numerical feedback models. 
Nevertheless, after trying many possible models and differ- 
ent parameters, they let us study very well the consequences 
of this assumptions on scales of interest. 
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6.3.1 Riddles 

A flaw we found in successful Bondi models is that they 
do not naturally generate X-ray cavities. Instead, the con- 
tinuous outflows carve a narrow and long tunnel along the 
z-axis. However, we have seen that continuous jets, with ram 
pressure almost equal to gas thermal pressure, are highly dis- 
rupted by turbulence, especially in central regions. Hence, 
the fragmentation of a feeble jet may produce generations 
of 'gentle' detached bubbles, which rise buoyantly, even in 
continuous feedback models. 

We propose that any kind of strong turbulence 
can indeed disrupt a m oderate power jet. For example, 
iMorsonv era!] (|2010') and' Pubois et all l|2010l 'l tried to start 
with a cluster atmosphere in non hydrostatic equilibrium, 
i.e. following a real cosmological evolution (local density fluc- 
tuations, merging, etc.), and they found a fragmentation of 
the AGN jets, helping the deposition of energy in the inner 
core. 

The creation of an unidirectional channel and subse- 
qu ent energy deposition at too large radii was pointed out 
bv lVernaleo &: Revnoldd (|200q ). However, if we perform long 
term evolutions, we have showed that the above mentioned 
turbulence and vorticity promotes mixing in the central ac- 
tive region, replenishing again the channel. Even in B0NDI2 
model, where a small narrow channel always stays open, in- 
stabilities and turbulence clearly heat the gas at the base of 
the jet, letting to cool only a moderate amount of gas in the 
near equatorial region. 



6.4 Best models confrontation 

It is interesting to flnd other important differences (or 
similarities) in successful models A and Bondi. The latter 
displays low velocities of the order 6 — 7 x 10^ km s~^, 
a value consistent w ith line- absorption observations (see 
ICrenshaw et al.l 120031 for a review, and other references in 
Section 1). Explosive models A tend to have also higher ve- 
locities, in a few events reaching ~ 10^ km s~^, more similar 
to a fast jet than an outflow wind. Both models have reason- 
able final injected energies, below a theoretical total energy 
of a BH with mass lO'', - 2 x lO'*^ erg. 

Apart from cavities and shocks (previously discussed), 
we produce other significant predictions, which are (or will 
be) comparable to X-ray observations. Iron abundance maps 
play a relevant role in tracing the outflow activity. The met- 
als, produced mainly by SNIa in the cD elliptical galaxy, are 
easily transported along jet-axis up to 150-200 kpc, creating 
an unmistakable asy mmetry, when the AGN is very active. 
The observations by iKirkpatrick et al.l (|2009l ) and Doria et 
al. (prep.) conflrm this behaviour. In the period of quies- 
cence turbulence and bulk motion will dominate the scene, 
smoothing and almost restoring the homogeneity. This could 
be another constraint in choosing between Bondi or not: in 
fact a continuous jet will tend to show often this marked 
asymmetry. 

Finally, another striking diversity is the clean depar- 
ture from hydrostatic equilibrium in models A, due to its 
more explosive nature. On the contrary, in Bondi type, mean 
fluctuations are very contained with errors below 8%. The 
consequence is, in the first case, a less precise mass deter- 



minati on using Eg. ([5l), as also found by oth er observational 
works l|Nagai et al.ll2007l : lRasia et al.lbood '). 

We conclude pointing out that the gaps of models A, 
can be replenished by the features of Bondi feedback, and 
viceversa (especially for cavities). Thus, an intriguing solu- 
tion may be a 'dual model', in which few energetic AGN 
outbursts (like models A3-4), perhaps triggered by accre- 
tion rates close to the Eddington limit, are superimposed (or 
alternated) to the weak activity induced by Bondi prescrip- 
tion. The first explosive mode will create large spheroidal 
bubbles, as those observed in real clusters, while the sec- 
ond quasi-steady outfiow will be the real sustaining pillar of 
the heating machine. This scenario would require a physical 
explanation of the alternation mechanism or why the two 
types of outfiow are diversified. 



7 CONCLUSIONS 

Overall we found that subrelativistic AGN outflows, pro- 
duced by two types of self-regulated feedback, are able to 
quench cooling flow for at least 7 Gyr and, at the same 
time, preserve the cool core appearance: 

(a) self-regulated feedback based on the instantaneous 
AA/cooi accreted, with mechanical efliciencies between 
~5 X 10"=* - IQ-^; 

(b) Bondi triggered feedback with e around 0.1, based on 
an almost continuous and steady outflow, generated by hot 
gas accretion. 

Both best models, (a) and (b), present primary merits 
concordant with X-ray observations: 

(la, b) cooling rate is reduced at least under 5% of the 
pure CF model; 

(2a, b) (mass-weighted) T{r) and n{r) profiles oscillate 
near the observed ones; 

(3a, b) total injected energy is always below 2 x 10®^ erg, 
under the limit of Ebh', 

(4a, b) mean Vjct'. from ~ 5 x 10'' to a few 10* km s~^. 

Intermittent (B) or continuous scheme with fixed ve- 
locity (C) are not consistent in the long term with some of 
these points, and hence rejected. 

In addition to these constraints we wanted to test fur- 
ther the best models. We warn that the following are supple- 
mentary predictions, but not main goals of our study: due to 
the need of a long term evolution, our simulations are lim- 
ited in resolution (~ kpc), and fine details could be altered 
or just missed by our numerical implementation. Marked 
differences between the two best models are: 

(5a) super-Eddington Mace: high power (~ lO**** erg s^^); 
(5b) sub-Eddington Mace: low power (~ 10*** erg s~^); 
(6a) cavities with high internal energy and shocked rims; 
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(6b) absence of a duty cycle and real inflated bubbles; 

(7a) asymmetrical transport of metals (SNIa) along 
jet-axis up to 200 kpc; subsequent gradient lengthened by 
turbulence; 

(7b) almost always asymmetrical iron maps; 

(8a) large deviation from hydrostatic equilibrium state 
(high turbulence); 

(8b) deviation from HE mass determination below 8% 
(moderate turbulence). 

Points (6) (caused by (5)) seem not to follow entirely 
observations (even if we need larger samples of clusters with 
AGN activity to reconstruct global history, and certainly 
higher resolved SB maps to study cavities). While the two 
models appear at first sight antithetical, it might be possible 
that they alternate during evolution, or that high power 
outbursts are just superimposed to a low Bondi feedback. 

Despite these riddles, purely mechanical AGN outflows 
promise to be good candidates, to uphold the wild and fre- 
netic 'dance' between heating and cooling. 
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